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r \ ■ ABSTRACT 

| To exploit synergies between the Herschel Space Observatory and next generation radio fa- 

. cilities, we have extended the semi-empirical extragalactic radio continuum simulation of 

Qh' Wilman et al. to the mid- and far-infrared. Here we describe the assignment of infrared spec- 

q \ tral energy distributions (SEDs) to the star-forming galaxies and active galactic nuclei, using 

$h ' Spitzer 24, 70 and 160/im and SCUBA 850/iin survey results as the main constraints. 

"^2 ' Star-forming galaxies dominate the source counts, and a model in which their far- 

d . infrared-radio correlation and infrared SED assignment procedure are invariant with redshift 

underpredicts the observed 24 and 70/im source counts. The 70,um deficit can be eliminated if 
the star-forming galaxies undergo stronger luminosity evolution than originally assumed for 
^ , the radio simulation, a requirement which may be partially ascribed to known non-linearity 

' in the far-infrared-radio correlation at low luminosity if it evolves with redshift. At 24/im, 

the shortfall is reduced if the star-forming galaxies develop SEDs with cooler dust and corre- 
spondingly stronger Polycyclic Aromatic Hydrocarbon (PAH) emission features with increas- 
ing redshift at a given far-infrared luminosity, but this trend may reverse at z > 1 in order 
, not to overproduce the sub-mm source counts. The resulting model compares favourably with 

■ recent BLASFresults and we have extended the simulation database to aid the interpretation of 

| Herschel surveys. Such comparisons may also facilitate further model refinement and revised 

predictions for the Square Kilometer Array and its precursors. 
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1 INTRODUCTION 

In Wilman et al. (2008) (hereafter W08) we presented a semi- 
empirical simulation of the extragalactic radio continuum sky pri- 
marily intended to aid the design of scientific programmes for the 
next generation of radio telescope facilities, culminating in the 
Square Kilometer Array (SKA). The simulation covers a sky area of 
20 x 20 deg 2 and contains radio-loud and radio-quiet active galac- 
tic nuclei (AGN) and star-forming galaxies out to redshift z = 20 
within a framework for their large-scale clustering. The full source 
catalogue - containing 320 million sources above 10 njy at five fre- 
quencies ranging from 151 MHz to 18 GHz - can be accessed at 
the SKADS Simulated Skies (S 3 ) database 1 under S 3 -SEX (semi- 
empirical extragalactic simulation). 

There are by necessity numerous uncertainties and limitations 
in the S 3 -SEX simulation, including but not limited to issues such 
as the form of the high-redshift evolution of the AGN and galaxies, 

1 http://s-cubed.physics.ox.ac.uk 



the lack of star-forming/ AGN hybrid galaxies, and the abundance 
of highly-obscured Compton-thick AGN. In so far as possible, flex- 
ibility was built into the simulations to allow post-processing to 
improve their accuracy as observations in the years ahead lead 
to improved constraints. Major advances in these areas are ex- 
pected from the far-infrared surveys to be performed by the Her- 
schel Space Observatory (Pilbratt 2008). To facilitate such com- 
parisons we have post-processed the S^-SEX simulation to cover 
these wavelengths. In this paper, we describe the recipes for as- 
signing infrared spectral energy distributions (SEDs) to the radio 
sources, using existing mid- and far-infrared results from Spitzer 
and sub-mm survey data from SCUBA as the primary constraints. 
We then present our predictions for Herschel surveys, bolstered 
by a comparison with results from the Balloon-bome Large Aper- 
ture Submillimetre Telescope (BLAST) which offer a foretaste of 
Herschel's capabilities. In keeping with the philosophy of the S 3 
project to maximise the degree of interaction between the user and 
the database, the infrared fluxes are also provided on the S 3 web- 
page from which users can generate simulation products for com- 
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parison with observations. The radio-infrared connection is of im- 
mense empirical and theoretical interest for extragalactic surveys, 
for the identification and follow-up of Herschel sources, and for 
probing the physics of the far-infrared-radio correlation and its 
possible evolution. The simulation can nevertheless also serve as 
a standalone Herschel simulation, for comparison with others such 
as the phenomenologically-inspired backward evolution models of 
Valiante et al. (2009) and Pearson & Khan (2009), and the model 
of Lacey et al. (2009) which is based on a semi-analytical galaxy 
formation model. 



2 THE ASSIGNMENT OF INFRARED SEDS TO THE 
RADIO SOURCE POPULATIONS 

The input radio source catalogue was obtained from the S 3 -SEX 
online database and negative evolution at high redshift was imposed 
using the 'default post-processing options' described in W08. Each 
radio source was assigned an infrared SED from various template 
libraries appropriate for star-forming galaxies and AGN, and the 
output flux density of each infrared model component for each 
galaxy [in log(F„(Jy)] was computed at observed wavelengths of 
24, 70, 100, 160, 250, 350, 450, 500, 850, and 1200 ^m. The wave- 
lengths of 24, 70 and 160 (im are those of the Spitzer MIPS in- 
strument, for which an abundance of extragalactic survey data is 
available to guide the construction of our model and to compare 
against its output (see the review of Soifer et al. 2008). Herschel 
will conduct surveys at 250, 350 and 500 /im with the SPIRE in- 
strument and at 70, 100 and 160/im with PACS. Longer- wavelength 
constraints are available from submillimetre surveys by SCUBA at 
450 and 850/im, and from the MAMBO instrument on the 30-metre 
IRAM telescope at 1200/im. The catalogued output flux densities 
are monochromatic values, with the exception of the Spitzer 24 and 
IQpm MIPS bands for which the modelled spectra were convolved 
with the MIPS bandpass transmission curves, following the Spitzer 
Synthetic Photometry Recipe 2 . This is due to the spectral complex- 
ity in the rest-frame 10/im region, particularly for the star formation 
component due to the presence of Polycyclic Aromatic Hydrocar- 
bon (PAH) features. 

The format of the output from the infrared extension, as it ap- 
pears in the S 3 -SEX online database, is described in Appendix A. 

2.1 Star-forming galaxies 

The population of star-forming galaxies in W08 comprises two 
populations, normal (or quiescent) galaxies and starbursts, identi- 
fied respectively with the low and high-luminosity Schechter func- 
tion components of the local 1 .4 GHz luminosity function of Yun, 
Reddy & Condon (2001). The entire population was evolved with 
pure luminosity evolution (1+z) out to z=1.5 (defined in an 
Einstein-de Sitter cosmology but adapted to the flat-A cosmology 
used for the simulation). The default post-processing option con- 
sists of negative space density evolution of the form (1 + z)~ 7 ' 9 
above z = 4.8. We stress that the terms 'normal galaxy' and 'star- 
burst' are merely convenient labels for these two components of the 
local luminosity function in our phenomenological model. Phys- 
ically speaking, the terms may not necessarily offer an accurate 
description of the nature of the star formation in these populations, 
especially beyond the local Universe. 

2 http://ssc.spitzer.caltech.edu/postbcd/cookbooks/synthetic_photometry.html 



The first step in assigning an infrared SED is to use the 
far-infrared-radio correlation for star-forming galaxies (eqn. 5 of 
Yun, Reddy & Condon 2001) in order to derive the rest-frame far- 
infrared luminosity, L(FIR), given the rest-frame 1.4 GHz luminos- 
ity, Li.4ghz- The relation is characterised by the "q" parameter: 

q = log[L(FIR, W)/3.75xl0 12 Hz]-log[L 1 .4GHz(WHz- 1 )],(l) 

for which we assume a gaussian distribution with /i = 2.34 
and a = 0.26 (Yun, Reddy & Condon). There is evidence that this 
relation holds out to redshift z > 1 (e.g. Appleton et al. 2004, Garn 
et al. 2009b), although such conclusions are sensitive to the as- 
sumed infrared SEDs and associated "k-corrections". Indeed, using 
recent BLAST stacking measurements, Ivison et al. (2009) reported 
evidence for a tentative decline in giR (defined as the logarithmic 
ratio of the the bolometric 8-1000/im infrared to monochromatic 
1.4 GHz radio fluxes), of the form q m oc (1 + z )-° ^±o.03 (gee 
also Kovacs et al. 2006). 

The far-infrared luminosity, L(FIR), is traditionally defined 
through a linear combination of IRAS 60 and 100/im flux den- 
sities (see Sanders & Mirabel 1996), which for an SED consist- 
ing of a superposition of modified blackbodies yields the 42.5- 
122.5/im luminosity to sub-percent accuracy (Helou et al. 1988). 
Given L(FIR), galaxies were assigned a model SED from the li- 
brary of star-forming galaxy templates assembled by Rieke et 
al. (2009). The latter consists of 14 SEDs uniformly spaced in total 
infrared luminosity (8-1000pm), L(TIR), from log L(TIR)( L Q ) = 
9.75 to 13.00. The templates were taken from the online version of 
Table 4 in Rieke et al. (2009) but they terminate on the short wave- 
length side at 5.26^im for those with log L(TIR) (L ) < 11 and 
at 4.02/im at higher luminosities. In order to predict 24pim fluxes 
for galaxies at redshift z > 3, we need to extend these templates to 
shorter wavelengths. For this purpose we employed the spectra of 
the individual local luminous (LIRGs) and ultraluminous infrared 
galaxies (ULIRGs) in Table 3 of Rieke et al. (2009), which extend 
down to 0.4^im. From these, we computed mean LIRG and ULIRG 
spectra and matched them onto the original Rieke templates below 
the cut-off wavelengths of 4.02 and 5.26pim. The ULIRG spectrum 
was used for the templates with log L(TIR)(L ) > 12, and the 
LIRG spectrum for the remaining templates. The resulting tem- 
plates should be considered as schematic below ~ 3/xm, as we did 
not attempt to model the stellar population and rest-frame optical 
extinction in a self-consistent fashion. 

The template SEDs were integrated over the appropriate wave- 
lengths ranges to yield the mapping between L(TIR) to L(FIR) 
shown in Fig. 1; the adopted relation between the two quantities 
is mildly non-linear: log L(FIR) = l.llog L(TIR) - 1.42. Based on 
the IRAS Bright Galaxy Sample (Soifer et al. 2007), Marcillac et 
al. (2006) assumed a linear relation: L(FIR) =1.91 L(TIR). In Fig. 2 
we show the Rieke templates scaled to a common L(FIR). 

At a given L(FIR), the local galaxy population exhibits a dis- 
tribution in 60-100/im colour or, equivalently, dust temperature, 
as characterised by Chapin et al. (2009). We use this distribution 
(Chapin et al. 2009, equations 8 and 9) to assign each star-forming 
galaxy a 60-100pim colour, C, defined as the ratio of the rest-frame 
monochromatic 60 and lOOpim flux densities: C = logLeo/Lioo- 
The template nearest in colour is then selected from the Rieke li- 
brary using the empirical C-L(TIR) conversion for these templates 
shown in Fig. 3, for which a functional fit is: 

C = 0.33tanh[l.l(logL(TIR) - 11.0)] - 0.23. (2) 

Having selected the template shape from the Rieke library 
by this process, the final step is to normalise it to actual value of 
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Figure 1. The symbols show the far-infrared luminosity (42.5- 
122.5 /im), L(FIR)(L ), of the star-forming galaxy templates of 
Rieke et al. (2009) as a function of total infrared luminosity, 
L(TIR). To assign SED templates to the simulated galaxies, we as- 
sume the non-linear relation given by the solid line: log L(FIR) = 
l.llog L(TIR) - 1.42. The dashed line shows the linear relation 
L(FIR) = 1.91 L(TIR) assumed by Marcillac et al. (2006). 
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Figure 2. The star-forming galaxy template SEDs scaled to a com- 
mon far-infrared luminosity using the L(TIR)-L(FIR) scaling in 
Fig. 1. With increasing L(TIR), the SEDs peak at shorter wave- 
lengths and exhibit weaker spectral features in the A < 20/im re- 
gion. 

L(FIR) originally specified by the far-infrared-radio relation. This 
is done using the non-linear L(FIR)-L(TIR) relation given in Fig. 1 . 
Chaplin et al. (2009) presented evidence for evolution in the colour- 
luminosity relation, but we initially assume that the local relation 
applies at all redshifts. 
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Figure 3. Colour of the Rieke templates [C = logL6o/Lioo] as a 
function of log L(TIR)(L ), with the fit given by equation (2) of 
section 2.1. 



findings of Hasinger (2008), who used a compilation of hard X- 
ray surveys to parameterise the fraction of obscured Compton-thin 
AGN (as a proportion of a total population comprising unobscured 
AGN and obscured Compton-thin AGN, but excluding Compton- 
thick obscured AGN) as a function of Lx (2-10 keV luminosity; 
ergs -1 ): 

h = [0.27 + /3(lo 9 L x - 43.75)]g{z), (3) 

where — —0.281 and g(z) is evolution of the form 
(1+z) 0,62 out to z — 2 and flat thereafter. The number of Compton- 
thick obscured AGN is not as well constrained, and in W08 we sim- 
ply boosted the space density of the Ueda et al. (2003) luminosity 
function by 50 per cent in a notional allowance for them, but left the 
issue open for subsequent refinement in post-processing. Ueda et al. 
showed that a reasonable match to the X-ray background results if 
the number of Compton-thick AGN simply equals the number of 
Compton-thin obscured AGN. We thus assume that the abundance 
of Compton-thick AGN population is a factor /cthick times that 
of the combined population of unobscured and Compton-thin ob- 
scured AGN, where /cthick = min[0.5,/2], with the upper bound 
of 0.50 hard-wired into the W08 simulation. After the removal of 
any excess Compton-thick AGN, the remaining sources in the W08 
catalogue are probabilistically identified as unobscured (type I), 
Compton-thin or Compton-thick obscured (type II) AGN. Finally, 
10 per cent of all sources are flagged for removal because W08 
included radio-loud AGN with a separate radio luminosity func- 
tion but this population is also implicitly present in the Ueda et 
al. luminosity function. The fractions of retained radio-quiet AGN 
which are unobscured, Compton-thin and Compton-thick obscured 
are shown as functions of Lx and redshift in Fig. 4. 



2.2 Radio-quiet AGN 

In W08 a population of radio-quiet AGN was incorporated using 
the Ueda et al. (2003) hard X-ray luminosity function and a re- 
lation between 2-10 keV and 1.4 GHz luminosities. To assign in- 
frared SEDs the first task is to split the population into subgroups of 
unobscured AGN, and Compton-thin and Compton-thick obscured 
AGN. This was not done in W08 but is now necessary because the 
infrared SEDs are sensitive to this classification. We start from the 



2.2.7 AGN emission 

For the retained sources, Lx is inferred from the radio luminos- 
ity using equation (2) of W08. Sources with Lx above and below a 
threshold of 3 x 10 44 erg s~ 1 are classified as quasars and Seyferts, 
respectively. Unobscured (type I) quasars are assigned at random 
one of the three templates for optically-selected QSOs compiled by 
Polletta et al. (2007) (termed QSOl, BQSOl and TQSOl, which 
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Figure 4. The fractions of retained radio-quiet AGN identified as unobscured AGN (solid lines), Compton-thin obscured (dashed lines) and 
Compton-thick obscured AGN (dotted lines), for intrinsic luminosity thresholds of Lx > 10 40 , 10 2 , 10 43 and 10 44 ergs _1 . Evaluated 
from a 4 deg 2 simulation area. 



differ in their relative optical/infrared flux ratios) 3 . The unobscured 
Seyferts are assigned SEDs from the online library of models com- 
puted by Nenkova et al. (2008) using the CLUMPY code for dust 
radiative transfer in a clumpy AGN torus. To narrow down the pa- 
rameter space of available models, we use the results of Thomp- 
son et al. (2009) which are based on analysis of Seyfert 1 AGN 
mid-infrared spectra. They prefer a range of CLUMPY parameter 
space defined by the parameters Y = 30, q = 1, a = 45 deg, 
m = 2, i = 30 deg, r v = 30 - 60 and A^o < 6 (see Nenkova et 
al. 2008 for the meaning of these parameters). 18 CLUMPY mod- 
els satisfy these constraints and were assigned at random to our 
type I Seyferts, with the SED consisting of the reprocessed torus 
emission and direct AGN broken power-law emission (defined as a 
power-law Fx oc A -1 ' 8 over 0.1-1/im with a Rayleigh- Jeans slope 
at longer wavelengths). 

For the type II Compton-thin AGN, we assigned one of the 
type 2 QSO templates given by Polletta et al. (2006) and Polletta et 
al. (2007) (termed QS02 and Torus, respectively) for the quasars, 
and one of the aforementioned CLUMPY models (but without the 
direct AGN emission component) to the Seyferts. 

For the Compton-thick AGN, we assigned half the popula- 
tion SEDs using the same prescription as for the type II Compton- 
thin sources; for the other half of the population, we reflect the 
growing body of evidence which suggests that a large proportion 
of Compton-thick AGN are obscured by an extended distribution 
of cold dust in the host galaxy and not by a nuclear torus (see e.g. 
Polletta et al., 2006, 2008; Martinez-Sansigre et al. 2006, Lacy et 
al. 2007). Accordingly we assign these sources an intrinsic type I 
or Compton-thin type II SED, extinguished by a dust screen with 
Ay = 4—27 mag (as found by Polletta et al. 2008 for obscured 

3 These SEDs were obtained from http://www.iasf- 
milano.inaf.it/~polletta/templates/swire_templates.html 



quasars). An extinction curve extending into the infrared is taken 
from McClure (2009). These AGN infrared SED components are 
normalised using the correlation between 12 micron and hard X-ray 
emission found by Gandhi et al. (2009), which appears to hold for 
obscured and unobscured AGN (including Compton-thick sources) 
and to extend into the quasar regime: 

logLMm = (-4.37 ± 3.08) + (1.106 ± 0.071)logL x , (4) 

where Lmir = A£a at 12.3/im in ergs -1 . Gaussian scatter 
of a = 0.36 is assumed on the relation. 



2.2.2 Star-formation 

There is abundant evidence for ongoing star formation in the host 
galaxies of radio-quiet AGN and this is likely to dominate the far- 
infrared SED. The consensus emerging from targeted observations 
of local Seyferts (Buchanan et al. 2006; Thompson et al. 2009) 
and studies of more distant samples and cosmological surveys (e.g. 
Silverman et al. 2009; Hernan-Caballero et al. 2009; Serjeant & 
Hatziminaoglou 2009) is that the star formation is generally higher 
in obscured AGN compared with their unobscured counterparts, 
and increases sub-linearly with AGN luminosity. We embody these 
findings with a star formation rate (SFR) oc V-^x(l + z) ' but 
assume different normalizations for the type I and II AGN. This 
luminosity and redshift dependence is similar to that measured 
for quasars by Serjeant & Hatziminaoglou (2009) and the lumi- 
nosity dependence matches that of local Seyfert Is (Thompson et 
al. 2009). 

For type I AGN, we normalise the relation using the mea- 
surements of type I SDSS AGN at z=0.15 by Kim et al. (2006), 
which have a median SFR of 0.5Meyr _1 . The median log 
L([OIII]A5007) of this sample, as inferred from Fig. 4(b) of Kim 
et al., is 41.4 (ergs" 1 ) which equates to log Lx = 42.6 using 
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the log Lx-log L([OIII]A5007) relation of eqn. (1) of Silverman et 
al. (2009). Thus for type I AGN: 

SFR = 0.63a/L x /10 43 (1 + z) 1 ' 6 M yr" 1 . (5) 

For the type II AGN (Compton-thin and Compton-thick), we 
base our results on the z = 0.1 SDSS type II AGN sample of 
Kauffmann et al. (2003), focussing on the sub-sample of strong- 
lined sources with log L([OIII])> 40.5 indicated by the red dots 
in Fig. 10 of Silverman et al. (2009). Referring once more to the 
log L([OIII]) distributions in Fig. 4(b) of Kim et al. (2006), we 
infer for this particular luminosity-limited sample a median log 
L([OIII])=41.1, which translates into log L x = 42.25 and hence 
boosts the pre-factor of our assumed SFR(Lx,z) relation from 0.63 
to 2M yr _1 . A scatter of 1 dex is assumed on the SFR(Lx,z) 
relation and the star-formation rates are assumed to decline as 
(1 + z)~ 7,9 at 2 > 4.8 following the assumed evolution of the 
global star formation rate density. SEDs for the star formation com- 
ponent are assigned from the Rieke et al. (2009) library, using the 
relation given by Yun, Reddy & Condon (2001) (their eqn. 13) to 
relate the star formation rate and 1.4 GHz radio luminosity; the far- 
infrared-radio relation (eqn. 1) was then used to assign L(FIR) and 
hence an appropriate template from the Rieke library using Fig. 1. 



2.3 Radio-loud AGN 

Radio-loud AGN were modelled by W08 using the Willott et 
al. (2001) 151 MHz luminosity function comprising a weakly- 
evolving low luminosity component and a rapidly-evolving high 
luminosity component, identified with FRI and FRII radio sources, 
respectively. We consider them separately for the purposes of in- 
frared SED assignment. 



2.3.1 FRIs 

There is scant evidence for the presence of intrinsically lu- 
minous optical-UV AGN in the radio population we asso- 
ciate with FRI sources (e.g. Vardoulaki et al. 2008). We there- 
fore model their infrared SEDs with a template representing 
host galaxy starlight only, which we produced using the GAL- 
SYNTH online interface to the GRASIL code developed by 
Silva et al. (1998). The models are based on their 'Ellipti- 
cal model', which gives the evolving SED (taking account of 
dust reprocessing) of a stellar population formed in a burst 
of 1 Gyr duration (see description of this particular model at 
http://adlibitum.oat.ts.astro.it/silva/grasil/modlib/modlib.html). We 
generated models with ages separated by 0.1 Gyr for 0.1-1.5 Gyr, 
and by 1 Gyr for ages 2-13 Gyr. We assume formation redshifts 
distributed uniformly in the range z= 10-20, and select the template 
nearest in age at the redshift of the galaxy under consideration. To 
normalise the SEDs, we draw upon the bivariate 1400 MHz-K- 
band luminosity function derived by Mauch & Sadler (2007) for a 
local sample: given a radio power, the absolute K-band magnitude 
is drawn from a Gaussian distribution of specified mean and vari- 
ance and this is used to normalise the SED. For simplicity, we as- 
sume the same (local) relation between K-band luminosity and ra- 
dio power at all redshifts. An artifact of this choice of model is that 
the mid-far infrared regions of the SED template increase sharply 
for ages below the assumed duration of the starburst (1 Gyr). This 
causes discontinuities in the mid- and far-infrared fluxes at z ~ 4 
where the galaxy age first drops below this figure. 



2.3.2 FRUs 

The FRIIs are first split into populations of quasars and radio galax- 
ies for which the nuclear emission is seen directly and through 
obscuration, respectively. Willott et al. (2000) found that for log 
L(151 MHz)[W/Hz/sr]> 26.5 the quasar fraction is 0.4, indepen- 
dent of redshift and luminosity, and around 0. 1 at lower radio lu- 
minosity. Taking into account the contamination by FRIs at lower 
luminosity, we simply assume that 40 per cent of our FRIIs are seen 
as quasars, i.e. those seen at viewing angle of < 53 degrees to the 
jet axis. To assign AGN and starburst infrared SEDs, we recall that 
quasars are, on average, more intrinsically luminous AGN than the 
radio galaxies (see e.g. Simpson 2003). To work with a parameter 
that is closely tied to the intrinsic strength of the nuclear emission, 
we convert from L(151 MHz) [W/Hz/sr] to L([OIII]A5007) [W] us- 
ing a relation from Grimes, Rawlings & Willott (2004) (GRW04): 

logL([OIII]A5007) = 7.53 + 1.0451ogL(151 MHz). (6) 

Eqn. (6) applies to the population as a whole, but following 
GRW04 we assume that the quasars and radio galaxies lie 0.3 dex 
systemtically above and below it, respectively, with an additional 
0.5 dex of scatter on each distribution. To assign an AGN com- 
ponent to the infrared SED, we draw upon the findings of Haas 
et al. (2008) who presented Spitzer 3-24^im photometry of 3CR 
steep-spectrum quasars and radio galaxies at redshift 1 < z < 2.5. 
The quasar mid-infrared SEDs are very uniform, roughly constant 
in vF v . There is a near proportionality between vF v {9> /^m, rest) 
and vF v (H9> MHz, rest) with 0.5 dex scatter. The radio galaxies are 
dominated by host galaxy starlight below 3/im (rest) and at longer 
wavelengths by a quasar nucleus reddened by up to Ay = 50 mag. 
The radio galaxies are a factor 3-10 less luminous in this spectral 
range than the quasars, which they attribute entirely to extinction. 
For a sample of 3CR sources at 0.4 < z < 1.2, Cleary et al. (2007) 
find that the quasars are 4 times as luminous as the radio galax- 
ies at 15/im, half of which they attribute to extinction, and half 
to synchrotron contamination in the quasars. However, their sam- 
ple includes a significant number of steep-spectrum compact and 
super-luminal quasars. They also assumed a spherically symmetric 
dust model, which is likely to be inadequate. 

We assume that the quasars and radio galaxies follow the same 
intrinsic relation between ^L„(8 pm rest) [W] and L([OIII]A5007) 
[W], derived from the quasar points in Fig. 1 of Haas et al. (2008) 
using eqn. (6) to convert from L(151 MHz) to L([OIII]A5007). We 
assume: 

logi/L„(8 ^m rest) = 3.3 + 0.95691ogL([OIII]A5007). (7) 

No intrinsic scatter is assumed on this, as the scatter in eqn. (6) 
can reproduce the scatter in Fig. 1 of Haas et al. (2008). Using 
this normalisation, we assign the quasars randomly to one of the 
three type I quasar templates from Polletta et al. (2007). The ra- 
dio galaxies are also assigned one of these template SEDs in an 
analogous manner, and then extinction of Av=0-40 mag is ap- 
plied. An elliptical host galaxy model from the GRASIL library 
is assigned as for the FRIs but in this case normalised using the 
K-z relation of Willott et al. (2003) with an additional offset of 
AK = -0.36(logL(151MHz)[W/Hz/sr] - 26.55) mag (follow- 
ing Willott et al. 2003, Fig. 3). 

In principle, the evolving elliptical template from the GRASIL 
library should self-consistently reproduce the far-infrared dust 
emission. To compare with this, we included in addition a far- 
infrared dust component to explicitly match the starburst properties 
to the sub-millimetre observations of radio galaxies and quasars 
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by Archibald et al. (2001), Willott et al. (2002) and Reuland et 
al. (2004). From Willott et al. (2002), the radio quasars appear to 
be a factor of several brighter than the radio galaxies at 850^im. 
In Grimes, Willott & Rawlings (2005), this was enshrined as a 
universal relationship between L([OIII]A5007) [W] and L(850^m) 
[W/Hz/sr]: 

logL(850Mm) = -2.5753 + 0.6821ogL([OIII]A5007), (8) 

with additional scatter drawn from a gaussian of a = 0.44. 
With this normalisation, we adopt a modified blackbody with /3 = 
1.95, T = 41 K, as found by Priddey & McMahon (2001) for high- 
redshift quasars. Archibald et al. (2001) found L(850/im (rest)) to 
be a strongly increasing function of redshift, but for our adopted 
SED parameters, the redshift dependence is a more modest (1 + 
z) 1 - 5 (see Reuland et al., Fig. 4); this dependence is added to 
eqn. (8) and normalised at z = 2.5. There is evidence for a de- 
cline in L(850/im (rest)) at z > 4, so in common with other parts 
of the simulation we assume a decline of the form (1 + z)~ 7 ' 9 
at z > 4.8. We note, however, that Rawlings et al. (2004) found 
any explicit redshift dependence to be small or absent. Although 
Willott et al. (2002) found tentative evidence for an anti-correlation 
between radio-source size and sub-mm luminosity for quasars, no 
such trend was found for the radio galaxy sample of Reuland et 
al. (2004) and we have chosen not to incorporate one. 

Finally, we remark that our treatment of the infrared SED is 
confined to reprocessed dust emission due to star formation and 
AGN heating. As a result, we have not endeavoured to incorporate 
an infrared synchrotron emission component in the beamed radio- 
loud AGN. If required, such a component could be added as an 
extrapolation of the radio SED given in the S database, following 
the beaming treatment described by W08. It would only have a sig- 
nificant effect on the predicted source counts at the brightest flux 
densities. 



3 COMPARISON WITH EXISTING OBSERVATIONS 

The prescriptions outlined in section 2 constitute our 'baseline 
model'. Here we compare its source counts and other derived quan- 
tities against observations. In Fig. 5 we show the normalised dif- 
ferential source counts at 24, 70 and 160/im for comparison with 
the most recent Spitzer MIPS survey results. To generate these 
plots we used a 4 deg 2 area of the input radio simulation, cou- 
pled with a shallow subset of the full 400 deg 2 radio simulation 
area to generate the bright flux ends (above 1.6 mjy at 24/im, 
30 inly at 70/im and 100 mJy at 160pim). The simulated counts 
fall decisively short of the observations below 1 mjy and lOmJy at 
24 and 70/im, respectively. There is a modest excess of simulated 
sources above ~ 100 mjy at 70 and 160pim, and above ~ 10 mJy at 
850/im (Fig. 6) where we compare with results from the SCUBA 
Half Degree Extragalactic Survey (SHADES; Coppin et al. 2006). 

Since the source counts are dominated by star-forming galax- 
ies, we show in Fig. 7 the quantity qm. = log(FmJ i^ocm) in 
the 24 and 70/im bands for the simulated normal galaxy popula- 
tion. Fir, and J^ocm are the simulated infrared and radio fluxes 
which would be observed in these bands, i.e. without the appli- 
cation of any K-corrections. The simulated sample is dominated 
by galaxies substantially fainter than those in current observa- 
tional samples, but our <j24 and qro values compare reasonably 
well with those measured locally (Appleton et al. 2004, Beswick 
et al. 2008 and Garn et al. 2009a,b). For a better comparison with 
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Figure 6. The thin continuous line shows the simulated integral 
source counts at 850/im for the baseline model using a 4 deg 2 re- 
gion of the input radio simulation; the contributions of normal and 
starburst galaxies are shown by the dashed and dot-dashed lines, re- 
spectively. The asterisks and thick lines show the measured counts 
and uncertainty region for the SCUBA Half-Degree Extragalactic 
Survey (Coppin et al. 2006). 

the observations, we applied cuts in flux density to yield the sub- 
samples also shown in Fig. 7: at 24/im, we selected galaxies with 
F24 > 0.3 mjy and F2o cm > 90 /iJy to mimic the selection of Ap- 
pleton et al. (2004), but we note that these authors only included 
sources with spectroscopic redshifts and appear not to discriminate 
between star-forming galaxies and AGN, biases which we cannot 
easily accommodate. At 70pim, we apply cuts of F70 > 6 mjy and 
i^ocm > 30 fily to match the selection of Seymour et al. (2009), 
who use photometric and spectroscopic redshifts and apply star- 
forming galaxy/AGN classification. The resulting sample shows 
only a very modest decline in 570 out to z = 2, consistent with 
the findings of Seymour et al. 



4 REVISIONS TO THE BASELINE MODEL 

We now adapt the baseline model to provide a better match to 
the existing infrared source count data. We focus entirely on the 
star-forming galaxies and perform the modifications in two steps, 
motivated by recent observational findings. Firstly, we use post- 
processing to modify the cosmological luminosity evolution of the 
star-forming galaxies to remove the shortfall in the 70pim source 
counts. Secondly, we allow the SED template assignment proce- 
dure to evolve towards cooler templates at higher redshift in an 
attempt to rectify the deficit in the 24/im counts. 

4.1 Modifications to the luminosity evolution of the 
star-forming galaxies 

In W08 we assumed that the radio luminosity function of the star- 
forming galaxy population undergoes pure luminosity evolution 
(PLE) of the form (1 + z) v with p=3.1 out to z=1.5, applied 
in an Einstein - de Sitter cosmology, but adapted to the flat A 
cosmology of the simulation. Since then, new observational con- 
straints on the evolutionary form have been published. At 70^im, 
Huynh et al. (2007) constrained the evolution out to z ~ 1 and 
found a degeneracy between luminosity and density evolution with 
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Figure 5. Source counts for the baseline model simulation (section 2) showing normal (asterisks) and starburst galaxies (cross), radio-quiet 
AGN (diamonds), and the total (plus signs) [FRIs and FRIIs are also included in the totals]. The large squares on the 24pm panel show 
the Spitzer counts of Papovich et al. (2004); the thick lines in the 70 and 160/im panels bound the observed counts in the S-COSMOS and 
SWIRE fields (Frayer et al. 2009). The plots were generated using a 4 deg 2 region of the W08 input radio simulation and the full 400 deg 2 
to generate the bright flux end (see section 2 for details). 



p — 2.78 in the PLE case; using several Spitzer surveys, Mag- 
nelli et al. (2009) found redshift evolution consistent with PLE with 
p = 3.6 ± 0.4 out to z ~ 1.3. Both these studies derived the evo- 
lution for a flat A cosmology, which for a given functional form of 
PLE, leads to a higher abundance of galaxies than present in the 
W08 simulation. The discrepancy is especially pronounced at and 
below the break in the luminosity function. We can, however, cor- 
rect for it by boosting the luminosity of the star-forming galaxies 
in the W08 catalogue by an amount AlogL(L, z). The boost re- 
quired to bring the W08 luminosity function into agreement with 
PLE of p — 3.1 to z = 1.5 in a A universe is shown in Fig. 8. [We 
note, however, that W08 assumed that the local luminosity func- 
tion is flat below 20.6 W/Hz, in which regime no luminosity boost 
can compensate for the deficit in space density; for continuity, the 
boost is accordingly held constant below this evolving break lumi- 
nosity]. After the application of this luminosity boost, the 24 and 
70/im source counts are as shown in Fig. 9. The 70 pm counts 
below lOmJy are now in satisfactory agreement with the observa- 
tional data but a deficit remains at 24/im, which we address in the 
next subsection. 



4.2 Evolution towards cooler dust templates and the final 
model 

There is growing evidence that the one-to-one correspondence be- 
tween infrared luminosity and template SED established in the lo- 
cal Universe breaks down at earlier cosmic epochs. To cite sev- 
eral examples, Symeonidis et al. (2009) showed that luminous and 
ultraluminous galaxies out to redshift z ~ 1 have much cooler 
dust distributions (i.e. SEDs peaking at longer wavelengths) than 
their local counterparts in the IRAS Bright Galaxy Sample (see 
also Symeonidis et al. 2008); and in the mid-infrared, Farrah et 
al. (2008) showed that ULIRGs at z ~ 1.7 have spectral features 
characteristic of starbursts with luminosities 10 10 ~ n Lq, rather 
than local ULIRGs. Such trends may reflect more spatially ex- 
tended star-forming regions at higher redshifts with lower dust op- 
tical depths. Based on a study of sub-mm galaxies at z ~ 2.5, 
Chapin et al. (2009) argued for luminosity evolution in the colour- 
luminosity relation of the form (1 + z) 3 , in the sense that the SED 
becomes cooler at higher redshift (see also Pope et al. 2006). 

Building on these findings, we now introduce phenomenolog- 
ical modifications to the SED template assignment procedure with 
the aim of reducing the remaining deficit in the 24/im source counts 
(Fig. 9). With reference to the Rieke templates in Fig. 2, we see that 
for a given L(FIR), evolution towards cooler dust templates with 
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Figure 9. 24, 70 and 160/im source counts for the modified model of section 4.1, in which the luminosity of the star-forming galaxies is 
boosted to mimic pure luminosity evolution in a A-universe. Symbols and observed data are as in Fig. 5. 



increasing redshift can significantly boost the observed 24/im flux 
relative to that at 70pim. Accordingly, the L(FIR) is scaled down by 
a factor f(z) for the purposes of selecting the far-infrared colour, 
C, from the C-L(FIR) distribution (see section 2.1). The cho- 
sen template is, however, still normalised to the original value of 
L(FIR). After some experimentation, we chose a function which 
evolves as f(z) = (1 + z) 2 ' 5 out to z = 1. Beyond z ~ 1, the 
evolution towards cooler templates must stop or go into reverse in 
order not to overproduce the 850^im counts. We chose a decline of 
the form (1 + z)~ 15 from 2 = 1 — 2, flat thereafter. Even with 
this in place, the starbursts still overproduce the bright end of the 
850/im source counts above lOmJy if we retain the 'default post- 
processing option', i.e. a (1 + z)~ 7 ' 9 decline in space density at 
z > 4.8. As shown in Fig. 10, one possible way to avoid this is 
to remove the entire starburst population (i.e. the high-luminosity 
component) at z > 1.5 and we adopt this solution for our final 
model. 

The removal of the high-redshift starburst component of the 
luminosity function implies that the evolving normal galaxy popu- 
lation in the simulation is sufficient to account for the star formation 
at these epochs. Although we used the terms 'normal galaxy' and 
'starburst' merely as labels to refer to the two Schechter function 
components of the local luminosity function, this trend may find 
some physical justification in the work of Obreschkow & Rawl- 
ings (2009). The latter authors showed that, with increasing red- 
shift, galaxy disks are smaller, leading to higher gas densities and 



pressures and hence higher molecular gas fractions and star forma- 
tion rates. In the local Universe, star formation rates in excess of 
~ 10 M© yr _1 can only occur in dense gas disks created by zero 
angular momentum gas in major mergers characteristic of starburst 
galaxies, e.g. the merger-triggered local ULIRGs. 

Even with the starburst cut-off in place, the modest (fac- 
tor ~ 2) excesses of simulated sources at the bright flux ends 
of all three Spitzer bands persist. In Fig. 12, we show the fiux- 
redshift planes for the normal galaxy and starburst populations for 
70 and 160/im fluxes exceeding 100 mjy; much of the excess nor- 
mal galaxy population occurs at redshifts z < 0.5, particularly at 
70pim. It may arise from a combination of effects, such as a de- 
parture from strictly power-law pure luminosity evolution or from 
an inaccurate treatment of the scatter on the far-infrared relation in 
eqn. (1) (i.e. from galaxies scattered into the high-L(FIR) gaussian 
wing of the "q" parameter distribution). In support of the first pos- 
sibility, Huynh et al. (2007) found no evidence for significant evo- 
lution at infrared luminosities below 10 11 Lq for redshifts z < 0.4. 

A modest (factor ~ 2) deficit also persists in the 24pim counts 
around 0.3 mjy, which suggests that the required template evolu- 
tion may be more extreme than assumed, possibly requiring a de- 
pendence on luminosity in addition to redshift. It is also possible 
that the Rieke templates underestimate the PAH strength in the 
higher redshift galaxies. 

Finally, we note that gravitational lensing biases (e.g. Negrello 
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Figure 7. The quantity qm — log{Fxa./ -RiOcm) is used to char- 
acterise the far-infrared:radio correlation for star-forming galaxies 
and is plotted here for the entire population of simulated normal 
galaxies (black dots) at 24 and 70/im, using the baseline model of 
section 2. The red crosses and diamonds are the sources remain- 
ing after the application of infrared and radio flux cuts to mimic 
the selection of Appleton et al. (2004) and Seymour et al. (2009), 
respectively. 

et al. 2007) are not accounted for in these simulations, but can be 
critical when the source counts are very steep. 

4.3 Comparison with observed redshift distributions 

Having achieved reasonable agreement with the Spitzer MIPS 
source counts, particularly at the faint end of the observed flux 
range, we now compare the redshift distributions. In Fig. 13, we 
show them for 24pm sources above flux limits of 0.08, 0.15 and 
0.3 mjy. The observational data were taken from Le Floc'h et 
al. (2009) and are based on photometric redshift analysis of 30,000 
sources in a ~ 1.68 deg 2 region of the COSMOS field. The sim- 
ulated distributions were derived from a 0.25 deg 2 region and are 
in good agreement with the observations. We note, however, that 
Le Floc'h et al. excluded galaxies brighter than Jab = 20 mag to 
minimise cosmic variance fluctuations at low redshift. As they also 
remarked, the small cosmic volume covered by the observations 
below z ~ 0.4 prevents any useful comparison in this regime. 

At 70jj,m we compare with results from the COSMOS and 
GOODS-N fields (Fig. 14) using a 4 deg 2 simulation area. The 
COSMOS results are based on photometric redshift analysis by 
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log L(1.4GHz) (W/Hz) 

Figure 8. The boost in luminosity, AlogL(L, z) (dex), which must 
be applied to the star-forming galaxies in the W08 catalogue to 
transform to pure luminosity evolution of the form (1 + z) 31 to 
z = 1.5 in a A-Universe. The lines show redshifts z=0.08 (solid), 
0.45 (short dashed), 0.83 (dotted), 1.2 (dot-dashed) and 1.5 (long 
dashed). 

Hickey et al. (in prep.) for a shallow field covering 1.43 deg 2 and 
a deep sub-region of 0.165 deg 2 . The deep field is complete to 
lOmly, but the shallow field is only 60 per cent complete at this 
flux level; in calculating dn/dz for the shallow field we compen- 
sated for this by weighting each source by l/c(/), where c(/) is 
the completeness level at its flux, as read off from Fig. 3 of Frayer 
et al. (2009). The shallow field dn/dz agrees better with the sim- 
ulation at z < 0.4, whilst at z > 0.75 the match to the deeper 
field is more favourable. Redshift distributions for 70/im COSMOS 
sources have also been compiled by Kartaltepe et al. (2009). 

At a deeper flux level of 2 mjy, we compare with the raw red- 
shift distribution derived by Huynh et al. (2007) for the 0.05 deg 2 
GOODS-N field, using a 0.25 deg 2 simulation area; there is again 
a deficit of observed sources but the GOODS-N source catalogue 
is only ~ 50 per cent complete at the flux limit (see Frayer et 
al. 2006a, Table 1), and the raw redshift distribution of Huynh et 
al. (2007; their Fig. 1) appears not to correct for this incomplete- 
ness. The effects of cosmic variance for the small field observed 
may also be important. 

At 850/im, the redshift distributions for flux limits of a few 
mjy (Fig. 10) have median redshifts at z = 2 — 3 and high-redshift 
tails extending beyond z = 4. Thus, even though the final model 
has eliminated the starburst (high-L) population at z > 1.5, it is 
nevertheless able to reproduce the population of sub-mm (SCUBA) 
galaxies which reside on average at 2 ~ 2.3, as deduced from spec- 
troscopic (Chapman et al. 2005) and radio-mm-FIR photometric 
measurements (Aretxaga et al. 2007). There is also observational 
evidence that the distribution extends to redshift z = 4 and beyond 
(Coppin et al. 2009). This reinforces our earlier statement that the 
terms 'normal galaxy' and 'starburst', which we use to label the two 
Schechter function components of the z = star-forming galaxy 
luminosity function, should not be invested with much physical sig- 
nificance beyond the local Universe. 



4.4 Possible evolution in the far-infrared-radio correlation? 

Our default assumption is that the luminosity boost introduced in 
section 4.1 applies equally to the radio and infrared luminosities 
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Figure 11. Spitzer MIPS source counts for the final model, after application of a luminosity boost (section 4.1) and template evolution (section 
4.2) to the star-forming galaxies and with the starburst population removed at z > 1.5; symbols and observed data are as in Fig. 5. 



and that it thus reflects a modification to the cosmological evolu- 
tion of the star-forming galaxy population as a whole. Indeed, as 
we show in Fig. 15, the 1.4 GHz source counts for the star-forming 
galaxies come into better agreement with the observations after the 
application of the luminosity boost. The relevant observations are 
taken from Seymour et al. (2008) and provide a break-down of 
the 1.4 GHz sources counts into contributions from AGN and star- 
forming galaxies, with the latter dominating below 0. 1 mjy. With 
reference to the radio source counts in Fig. 4 of W08, the increase 
in the star-forming galaxy contribution brings the total radio source 
count at 1.4 GHz into better agreement with the observations at 
this flux level. In light of this, the luminosity boost should be ap- 
plied when using the simulated radio catalogues, as described in 
Appendix A. 

Nevertheless, we cannot exclude the possibility that some con- 
tribution to AlogL(L, z) may arise from redshift evolution and 
non-linearity in the far-infrared:radio correlation. Observational ev- 
idence for redshift evolution is at present limited, as discussed in 
section 3, but strong upward evolution in the relation is expected at 
z > 3 due to increased inverse Compton scattering off the cosmic 
microwave background (e.g. Murphy 2009). In the local Universe, 
there is evidence that the relation may become non-linear at low 
luminosities (e.g. Condon 1992; Yun, Reddy & Condon 1992; Best 
et al. 2005). Best et al. (2005) found that below log L(1.4 GHz) 
= 22.5 W Hz -1 the measured local radio luminosity function falls 
below that derived in the far-infrared. Such non-linearity could par- 



tially account for the luminosity boosting, but it would need to 
evolve strongly with redshift in order to match that required in 
Fig. 8. 

At higher redshift, major contributions in this area are ex- 
pected from Herschel surveys, as suggested by initial results from 
BLAST (Ivison et al. 2009). To compare with the latter we show in 
Fig. 16 the simulated quantity gm,, defined analogously to eqn. (1) 
but with the far-infrared luminosity replaced by the rest-frame 8- 
1000/im luminosity, L(TIR), (eqn. 4 of Ivison et al). The simulated 
values are in good accord with the value of qm, — 2.41±0.2 for the 
250/im-selected sample of Ivison et al. (2009), with no evidence of 
redshift dependence. Ivison et al. also determined qm for a sam- 
ple derived by stacking at the positions of 24/im sources and found 
tentative evidence for a decline of the form (1 + 2 )-° 15 ± 03 5 but 
we do not attempt to replicate the selection effects associated with 
the definition of this particular sample. The lack of redshift depen- 
dence in the simulated qm primarily reflects the redshift-invariance 
we assumed for the far-infrared-radio relation in eqn. (1), in con- 
junction with the only mildly non-linear dependence of L(TIR) on 
L(FIR)inFig. 1. 

For completeness, we also show in Fig. 16 the quantities q24 
and qro as functions of redshift, for comparison with those of the 
baseline model shown in Fig. 7. Reflecting the imposed template 
evolution, <j24 increases more strongly with redshift out to z = 1 
than in the baseline model, whilst 970 remains effectively constant 
out to 2 = 1. 
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Figure 13. Comparison of the predicted 24^im redshift distributions for the final model of section 4.2 (solid lines) in the central 0.25 deg 2 
of the simulation with those measured in a 1.68 deg 2 region of the COSMOS field by Le Floc'h et al. (2009). The latter authors excluded 
galaxies brighter than Jab = 20 mag to minimize cosmic variance fluctuations at z < 0.4. 



5 PREDICTIONS FOR HERSCHEL SURVEYS 



We now present some predictions for forthcoming Herschel sur- 
veys and a comparison with results from the BLAST mission which 
flew a prototype of the SPIRE instrument on a telescope half Her- 
schel's size. In Fig. 17, we show predicted integral source counts 
for surveys at 100, 250, 350 and 500/mi which correspond closely 
with those predicted by the backward-evolution model of Valiante 
et al. (2009). Also shown in this figure are normalised differential 
source counts at 250, 350 and 500/im, which compare reasonably 
well with the BLAST counts as determined from the 'P(D)' analy- 
sis of Devlin et al. (2009) (see also Patanchon et al. 2009). At 250 
and 350/im, the model marginally exceeds the BLAST measure- 
ments in the highest flux bins, as also hinted at in the 160/im com- 
parison above ~ 100 mjy. In Fig. 18 we compare with the 'com- 
plete' redshift distribution derived by Dunlop et al. (2009) for the 
150 arcmin 2 GOODS-S field at 250/im; the observed sample con- 
sists of 20 sources with redshifts down to a nominal flux limit of 
35 mjy, although small number statistics and cosmic variance limit 
the scope of the comparison. 



schel Key Programme 4 surveys: (i) the ultra-deep 0.012 deg 2 'pen- 
cil beam' GOODS-Herschel field 5 ; (ii) Levels 1 and 5 of the Her- 
schel Multi-tiered Extragalactic Survey (HERMES 6 ); (iii) the wide 
(~ 600 deg 2 ) but shallow ATLAS survey (Eales et al. 2009). Ta- 
ble 1 shows their area coverages, multi-band sensitivities, predicted 
numbers of galaxies and the fractions at z > 1 and z > 2. Our 
dn/dz predictions were derived from a 0.25 deg 2 simulation area 
and scaled to the required survey areas, except for the ATLAS and 
HERMES Level 5 predictions for which we used a 4 deg area . The 
predicted redshift distributions correspond quite well with those of 
Lacey et al. (2009) (L09), despite the differences in our simula- 
tion methodologies. For GOODS-Herschel, our predicted redshift 
distributions exhibit broader peaks than those of L09, extending be- 
yond 2 = 1; similarly, for HERMES Level 1 and 5 at 250/im our 
distribution peaks just above z = 1 and that of L09 just below. At 
100 and 160/im in HERMES Level 5, we predict broad redshift dis- 
tribution with a tail extending beyond z = 1 whereas L09 predict 
distributions strongly peaked at z — 0.2. For the ATLAS survey, 
the distributions of L09 exhibit marginally higher peaks than ours 
at 2 < 0.2; the L09 distributions cut-off sharply above z ~ 1.5, as 



Although the S 3 online interface allows the user to query 

the database to generate redshift distributions for arbitrary ra- 4 http://herschel.esac.esa.int/Key_Programmes.shtml 

dio and infrared flux cuts and survey areas, we follow Lacey et 5 http://herschel.esac.esa.int/Docs/KPOT 

al. (2009) and show in Fig. 16 dn/dz plots for a selection of Her- 6 http://astronomy.sussex.ac.uk/~sjo/Hermes 
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Table 1. A selection of Herschel Key Programme surveys for which we provide dn/dz plots in section 5. 



Survey 


Area 
(deg 2 ) 


Band 
(fim) 


Depth 
(mJy) 


Total* 


f(z>l)t 


f(z>2)ft 




U.U1Z 




ft 
U.D 




U.Ol 


14 


Ultra-deep 




160 


Q 


460 


69 


25 


Hermes Level 1 


0.11 


250 


4.2 


1850 


0.69 


0.24 






350 


5.7 


1080 


0.75 


0.29 






500 


4.9 


680 


0.82 


0.37 


Hermes Level 5 


27 


100 


27 


11440 


0.24 


0.0018 






160 


39 


13420 


0.41 


0.033 






250 


14 


85700 


0.63 


0.18 






350 


19 


27000 


0.69 


0.24 






500 


16 


10470 


0.79 


0.34 


ATLAS 


600 


100 


67 


59700 


0.15 


1.6E-04 






160 


94 


57170 


0.30 


8.3E-03 






250 


46 


165100 


0.51 


0.10 






350 


62 


28530 


0.51 


0.12 






500 


53 


5720 


0.63 


0.19 



* Total number of simulated galaxies in survey area at z = 1 — 4 
f Fraction at z > 1 
ft Fraction at z > 2 



do ours, due to the step function cut-off in the starburst population 
above z = 1.5. 



able model parameter space. Data products from our simulation 
are available on the S 3 website 7 and we expect them to serve as a 
valuable resource for the interpretation of Herschel surveys. 



6 CONCLUSIONS 

We have post-processed the S 3 -SEX semi-empirical simulation of 
the extragalactic radio continuum sky (W08) to make predictions 
for Herschel surveys at far-infrared wavelengths. Existing obser- 
vations in the mid-infrared with Spitzer at 24, 70 and 160/im and 
at 850/im with SCUBA, together impose strong constraints on the 
assignment of infrared SEDs to the star-forming galaxies as a func- 
tion of redshift. Our principal findings, incorporated into the final 
model, are as follows: 

(i) In order to match the 70/itn counts, the star-forming galax- 
ies are required to undergo stronger luminosity evolution than as- 
sumed by W08 for the radio simulation. It may be that W08 sim- 
ply used an inaccurate evolutionary prescription based on the in- 
formation available at that time; alternatively, it could be that there 
is genuine differential evolution between the far-infrared and ra- 
dio populations, as a result of an evolving non-linearity in the far- 
infrared:radio correlation which is already apparent locally at low 
luminosity . Results from Herschel and SKA precursors will resolve 
this issue. 

(ii) From the local Universe to redshift 2 = 1, star-forming 
galaxies are required to develop progressively cooler SED tem- 
plates at fixed L(FIR) once the latter has been set by the far- 
infrared-radio correlation; the rest-frame 60-100/im colour is as- 
sumed to evolve as logLeo/Lioo ~ (1 + z) -2 ' 5 ; beyond z = 1, 
this evolution must go into reverse in order not to overproduce the 
850/im source counts. 

Our chosen model compares favourably with recent far- 
infrared survey results from BLAST. This inspires confidence in 
using it to make predictions for Herschel Key Programme surveys. 
Our predicted source counts and redshift distributions correspond 
closely with those of Valiante et al. (2009) and Lacey et al. (2009), 
respectively, despite the clear differences in the methodologies of 
our models. This suggests that the combination of Spitzer mid- 
infrared and 850/itn data impose strong constraints on the allow- 
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Figure 10. Upper panel: 850pim source counts for the final model 
of section 4.2 as determined from a 4 deg 2 simulation area; the 
observed SHADES counts are the same as in Fig. 6; the dashed line 
shows the contribution of the normal galaxies and the dot-dashed 
line represents the starbursts, assuming a sharp cut-off in the latter 
population at z > 1.5; the dotted line shows the contribution of 
the starbursts under the assumption of the 'default post-processing' 
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panel: 850/im redshift distributions for flux limits of 2.4 and 5 mJy, 
assuming the z > 1.5 cut-off in the starburst population. 
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Figure 14. Comparison of the simulated redshift distribution for the 
final model of section 4.2 (solid lines) with observations at 70/im. 
Upper panel: results from the COSMOS field from Hickey et al. (in 
prep) for S(70/im)> 10 mJy (dashed and dotted histograms are for 
the deep and shallow fields, respectively - see section 4.3 for de- 
tails); lower panel: the observed data (histogram) are the raw red- 
shift distribution of Huynh et al. (2007, Fig. 1) for the GOODS-N 
field above 2 mjy, for which the 70pim source catalogue is, how- 
ever, only about 50 per cent complete. 
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Figure 15. Differential 1 .4 GHz source counts for the star-forming 
galaxies in the simulation with (dashed line) and without (continu- 
ous line) the application of the luminosity boost introduced in sec- 
tion 4.1. The observational data points are the star-forming galaxy 
source counts measured by Seymour et al. (2008). 



16 columns of which contain the existing output from the W08 ra- 
dio simulation. Each row of this table refers to an individual galaxy; 
multi-component radio sources (i.e. FRI and FRII sources) are de- 
noted by their core position and integrated radio fluxes. A descrip- 
tion of the supplementary columns is given in Table Al. 

In addition to the mid- and far-infrared flux densities (24- 
1200/im), we also provide 2.2/im K-band magnitudes. For the 
radio-loud AGN, the K-band magnitude is used to normalize the 
SED; for the radio-quiet AGN and star-forming galaxies, the ri- 
band magnitude should be considered as schematic as no attempt 
was made to model accurately the stellar population and dust ex- 
tinction at rest-frame optical wavelengths, even though the SED 
templates extend into this regime. 

The RQ-AGN classification flag applies only to the RQ- 
AGN and specifies whether the source is classified as unobscured, 
Compton-thin or Compton-thick obscured, or whether it is part of 
the 'excess' population which is filtered out (as described in section 
2.2). 

The Space-density filter flag indicates which sources have 
been filtered out due to the imposed cut-off in the space density 
at high-redshift. The infrared fluxes of all filtered-out sources are 
set to -999. 

Alog(L,z) is the luminosity boost which has been applied to 
the star-forming galaxies to generate the infrared fluxes, over and 
above the level set by the far-infrared:radio correlation. Note that 
the tabulated radio fluxes for the star-forming galaxies are the origi- 
nal values from W08, even though the boost may also be applicable 
to them, as discussed in sections 4. 1 and 4.4. 



APPENDIX A: DESCRIPTION OF DATABASE FORMAT 

The simulated infrared fluxes are available from our online inter- 
active database at http://s-cubed.physics.ox.ac.uk/s3_sex. They are 
provided as supplementary columns to the Galaxies Table, the first 
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Figure 16. Upper panel: the bolometric (8-lOOO^tm) gm for the sim- 
ulated star-forming galaxies. The dashed lines represent the ±3cr 
range measured by Ivison et al. (2009) for a 250/im sample. The 
continuous line shows the tentative evolution measured by Ivison 
et al. from stacking measurements. Middle and lower panels: 524 
and §70 for the final model, with symbols defined as in Fig. 7. 



Figure 17. Upper panel: predicted Herschel integral source counts at 
100 (black solid line), 250 (red dotted line), 350 (green dashed) and 
500/iin (blue dot-dashed line); lower panel: in differential format 
for comparison with the BLAST counts from Table S2 of Devlin et 
al. (2009), which are shown by the points with error bars (asterisks, 
250/im; squares, 350/im; triangles, 500/xm). 



BLAST GOODS-S (250 micron) 




Figure 18. A comparison of the redshift distribution of BLAST 
250/iin sources brighter than 35 mjy in GOODS-S (Dunlop et 
al. 2009) (histogram) with the simulated distributions for flux limits 
of 25, 30, 35, 40 and 50 mjy (solid lines, top to bottom), determined 
from a 0.25 deg 2 simulation area. 
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Table Al. Supplementary columns in the S 3 -SEX Galaxies Table derived 
from the infrared extension to the radio simulation. 



Column 


Attribute description 


[1-16 


Existing radio simulation output] 


17 


logi(j(24^/m flux density) [Jy] 


18 


logio(70^im flux density) [Jy] 


19 


lr»(Ti /i ( 1 00/ *m flnv Hpncitv*i tTvl 
\ M 11UA V J L YJ 


20 


logio(160/im flux density) [Jy] 


21 


logi o (250^im flux density) [Jy] 


22 


logio(350/tm flux density) [Jy] 


23 


logio(450/jm flux density) [Jy] 


24 


logio(500/im flux density) [Jy] 


25 


login (850/jm flux density) [Jy] 


26 


logio(1200/^m flux density) [Jy] 


27 


2.2/rni K-band magnitude 


28 


RQ-AGN classification flag * : 




1=TYPE 1 (unobscured), 




2=TYPE 2 (Compton-thin obscured), 




3=TYPE 2 (Compton thick obscured) , 




-l=Source excluded 


29 


Space-density filter flag: 




(-1= filtered out; 0=retained) 


30 


Alog(L,z) f 



* Applies to radio-quiet AGN only, 
f Applies to star-forming galaxies only. 
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Figure 19. Predicted redshift distributions for a selection of Her- 
schel Key Programme surveys (listed in Table 1) at 100/im (black 
solid lines), 160^im (purple dot-dot-dashed lines), 250/im (red dot- 
ted lines), 350/^m (green dashed lines) and 500/im (blue dot-dashed 
lines) 



